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Abstract 

In this work we present a method, based on the use of Bernstein poly- 
nomials, for the numerical resolution of some boundary values problems. 
The computations have not need of particular approximations of deriva- 
tives, such as finite differences, or particular techniques, such as finite 
elements. Also, the method doesn't require the use of matrices, as in res- 
olution of linear algebraic systems, nor the use of like-Newton algorithms, 
as in resolution of non linear sets of equations. An initial equation is 
resolved only once, then the method is based on iterated evaluations of 
appropriate polynomials. 

1 Basic concepts 

Let B n> i the i-ih n-degree Bernstein polynomial (see pQ): 

s»,i(*) = (?)t i (i-*) (n - i) , te[o,i] (l) 

Let Pi — (pi, qi) £ R 2 , with < i < n. Then we can consider, from [0, 1] to R 2 , 
the Bezier curve (see [2]) spawned by the array of points (Pi)o<i<«: 

n 

bz(t) = Y J B n At)V l (2) 

j=0 

Let y — f(x) a real differentiable function define on an interval [a, b) of K, and 
let (4>(t), (j)(t)), with t £ [0, 1], a possible parametric representation for /. Then, 
if t £ [0,1], we have x = tp{t), y{x) = f(ip{t)) — (f>(t) and, if x = t/^o), 
4>'{h) = f(x W{t ) (see El). Therefore, if V'(*o) ¥= 0, 



i 



If / € C k ([a, b\), than the Bezier curve 



bz{t) = Y,B n M( l -J{^]) (4) 



and its derivatives until the A:-th order are uniformely convergent to / and its 
derivatives for n — > +00. 

For our purposes, we'll transform a differential equation _F(?/ fc ) (x), x) = 
x G [a, 6], from cartesian to parametric representation. As example, using the 
above notations and relations, the first degree ODE y'(x) = x + y 2 (x), supposed 
tp'it) ^ in (0, 1), can be written in this form: 

ffl = m + m)f (5) 

Boundary values y(a) — y a , y(b) = yb can be written in the form 0(0) — 4> a = Vai 
0(1) = (j)(b) = yb, where a = ip(Q) and b = ip(l). 

A result on the theorical resolution of an initial value problem for a particular 
ordinary differential equation is given in (@|). 



2 The algorithm of Bernstein polynomials 

Suppose we want a numerical solution for a BVP writable in the form 

F(<t>W(t),ii>(t))=4>(t), te[o,i], 0(0) = 00,0(1) (6) 

where < k < N for some N G N+. Let A the array of the control points 
{(a,ya), (b,yb)}, and dt the grid step for the parameter t. 

1. As first part of the algorithm, let (p,q) a pair of unknowns, and let 
621 = (ijji(t), 0i (i)) the Bezier curve spawned by the control points array A = 
{(a, y a ), {p, q), (b, yb)}- Then we compute the values of p and q using the follow- 
ing two equations, obtained from the differential equation and from the bound- 
ary conditions respectively at t — and t=l: 



F(^>(0),a)=y a 
F(4 k \l),b) = y b 

This system can have more then one real solution; in this case we choose the 
pair with the value of p next to 0.5(6 + a), the medium point of the cartesian 
interval [a, b\. We'll call the pair (p, q) the pivot point. If the system have not 
real solutions, this first step has no output. 

The pivot is useful to obtain a first raw graphic profile of the solution y{x), 
since the line from (a, y a ) to (p, q) is tangent to the Bezier curve at (a, y a ), and 
similarly for (b, yb) and (p, q) (see 0). 



2 



2. As second part of the algorithm, we open a loop starting from m = 1. We 
evaluate the last Bezier curve bz m — (ip(t), 4>{t)), spawned by the actual control 
points array A, at the values to = (mdt) and ti = (1 — mdt) of the parameter t. 
Note that tp(t) and 0(t) are polynomials. In this way we obtain the two points 
(^(to), (f>{to)) and (ip(ti), (f>(ti)). Then we can compute the local errors 



eo=F(^*>(t o ),V(«o))-0(to), ei = F(0W(ti),^(ti))-0(ti) (7) 

These errors are an estimate about the approximation of the Bezier curve bz m 
as solution of the BVP @. With the values eo and ei, we can use a predefined 
criterion for deciding to exit from the loop or not. In the present work we use 
the following simple but effective criterion: let s new = 0.5(ei + e-z) and s u the 
corresponding value from the previous loop step. If |s„ eu ,| > \s id\, then we exit 
from loop, else we insert the new points (ip(to), 4>(to) + eo),(^>(ii), </>(£ 1) + ei) in 
the correct place of the control points array A, and perform a new step. 
In any case, the loop finishes when a predefined maximum number of iterations 
is reached or when the product mdt is greater then the half of the interval [0, 1]. 

This technique can be considered a particular case of the collocation method (sec 
[S]) for the numerical solution of differential problems. 



3 A BVP with a linear differential equation 

As first application of the algorithm, we consider the following BVP: 

y'{x)=y{x) Vxe(0,2), y(0) = 1, y(2) = e 2 (8) 

which has the analytical solution y — e x . According to J2J), the problem has the 
following form when referred to the new variable t S [0, 1]: 

|tt| = <K*) Vie (0,1), #0) = 1, q>(l) = e 2 (9) 

The first step of the algorithm gives as pivot the point (0.56344, 1.68501). With 
a parametric step equal to 0.1, the algorithm exits only after 3 iterations, with a 
set of 9 control points. Then the analitical approximated solution of the problem 
is a Bezier curve of 8-th degree. In the following table we compare the numer- 
ical values between the approximated solution and the exact solution e x ^ for 
t e {0.1i} <i<io- 
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CXIJ LJ± UA 




0.0 


1.0 


1.0 


0.1 


1.17821 


1.19321 


0.2 


1.38942 


1.42075 


0.3 


1.75235 


1.77812 


0.4 


2.35041 


2.34651 


0.5 


3.18769 


3.15902 


0.6 


4.17184 


4.14995 


0.7 


5.14898 


5.14912 


0.8 


5.98714 


5.99204 


0.9 


6.67677 


6.67179 


1.0 


7.38906 


7.38906 




0.5 1 1.5 2 

Figure 1 : The continuous line is the graphic of the exact solution, while dots are the 
cartesian pairs of the numerical solution as in previous table. 

In Fig.l we compare the graphics of the two solutions. 

4 A BVP with a non linear differential equation 

In this section we consider the following BVP: 

-\(v\x)) 2 +A = y{x) Vxe(-1,3), y(-l) = 3, 2/(3) = -5 (10) 

which has the analytical solution y = —x 2 + 4. In this case the solution has a 
local maximum in the interior of the domain. The problem has the following 
form when referred to the new variable t G [0, 1]: 
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The first part of the algorithm gives as pivot the point (1.0, 7.0). With a para- 
metric step equal to 0.1, the algorithm exits only after 3 iterations, with a set of 
9 control points. In the following table we compare the numerical values between 
the approximated solution and the exact solution — x 2 (i)+4 for t £ {O.H}q<;<io. 



t 


approx 


exact 


-1.0 


3.0 


3.0 


-0.721242 


3.43889 


3.47981 


-0.439246 


3.80217 


3.80706 


-0.0632671 


4.02617 


3.996 


0.429831 


3.8453 


3.81525 


1.0 


3.02132 


3.0 


1.57017 


1.56462 


1.53457 


2.06327 


-0.226895 


-0.257071 


2.43925 


-1.95482 


-1.94992 


2.72124 


-3.44608 


-3.40516 


3.0 


-5.0 


-5.0 




Figure 2: The continuous line is the graphic of the exact solution, while dots are the 
cartesian pairs of the numerical solution as in previous table. 

In Fig. 2 we compare the graphics of the two solutions. 

5 A more subtle BVP 

In this section we consider the following BVP: 

y'(x) = x + y 2 (x) Vxe (0,0.9), y(0) = 1, 3/(0.9) = 32.725 (12) 

which cannot be solved in terms of elementary functions (see [H]). The built-in 
function NDSolve of the package Mathematica 4-0 (see [7]) is not able to solve it 
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directly, but we have numerically resolved the corresponding initial value prob- 
lem for the condition y(0) — 1, and then found the value for y(0.9). Note that, 
in the interval [0, 0.9] considered, the derivative y' is positive, hence the solution 
y is increasing. With a parametric step equal to 0.15, the algorithm of Bernstein 
polynomials exits after 2 iterations, with a set of 7 control points. In the follow- 
ing table we compare the numerical values between the approximated solution 
computed with this method (Bps) and the approximated solution computed 
using NDSolve. 



X 


Bps 


NDSolve 


0.0 


1.0 


1.0 


0.13872 


1.87331 


1.17177 


0.276777 


2.77787 


1.43186 


0.418191 


3.77639 


1.85729 


0.554867 


4.80733 


2.57996 


0.674855 


5.94924 


3.84768 


0.768333 


7.55906 


6.12075 


0.831326 


10.2845 


10.0499 


0.867159 


14.9504 


15.7566 


0.885645 


22.3191 


22.256 


0.9 


32.725 


32.725 



We note that locally the differences of values between the two methods can be 
large, but globally our method gives the same qualitative informations as the 
Mathematica internal function, in particular the growth in the interval [0, 0.9] 
and the presence of a value xq > 0.9 such that lim^a;,, y(x) = +oo. In Fig. 3 
we compare the graphics of the two solutions. 




0.2 0.4 0.6 0.8 



Figure 3: The continuous line is the graphic obtained by the Mathematica function 
NDSolve, while dots are the cartesian pairs of the solution provided by the method of 
Bernstein polynomials. 
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6 Technical notes 



All the computations have been made using the software Mathematical ver. 4.0, 
of Wolfram Research, on an Intel Xeon 3.2 GHz system. For an introduction 
to the use of this package in symbolic and numerical resolution of differential 
problems, see (0). 

The basic Bernstein polynomials are build with the instructions 

Bernstein [i_,n_] : =Binomial [n, i] t~i (1-t) " (n-i) ; 
Bernstein[nJ : =Table [Bernstein [i ,n] ,{i,0,n}] 

The symbolic derivative of Bernstein polynomials are computed with 

D [Bernstein [i ,n] ,t] 

The Bezier curve are obtained using dot product: 

Bezier [controlPoints_] : =Bernstein [Length [controlPoints] -1] . controlPoints 

The computation of the pivot point is made using the built-in function NSolve 
for the resolution of a system of two equations. The evaluations of the Bezier 
curves at t- values are made using the transformation rules / . operator. 

Electronic copies of the notebooks used for the numerical examples described in 
this work can be requested at the mail address gianluca.argentini@gmail.com. 
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